#!/bin/env python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from scipy.stats import norm

tfkinds = {1:'TC+SF',2:'TC',3:'SF',4:'Ohters'}

regname={
   1: "E.Japan",
   2: "C.Japan",
   3: "W.Japan",
}

sig_level=0.05
fit_norm = True # normal distribution

if __name__ == "__main__":

#   #---Observations (HGRID2_Debus_obsonly3, i063)
    df_obs = pd.read_json("./Data/anomalousdays_Observations.json")


#   #---HiFLOR AllForc each ens 1977-2015
    df_all = pd.read_json("./Data/anomalousdays_HiFLOR_AllForc.json")
    print (df_all)

#   #---HiFLOR NatForc each ens 1977-2015
    df_nat = pd.read_json("./Data/anomalousdays_HiFLOR_NatForc.json")
    print (df_nat)

    tkind = "ALL"
    xranges = {
        "ALL":[-0.3,0.5],
    }

   #--plot
    fig, ax = plt.subplots(figsize=(8.0,5.0))
    
    df_all[df_all["Kind"]==tkind].plot.hist(density=True, ax=ax, color="red", alpha=0.3, label="AllForc")
    if fit_norm: 
        sns.distplot(df_all[df_all["Kind"]==tkind]["Trend"].to_numpy(), ax=ax, fit_kws={"color":"red","linewidth":3}, fit=stats.norm, hist=False, kde=False)
    else:
        sns.distplot(df_all[df_all["Kind"]==tkind]["Trend"].to_numpy(), ax=ax, kde_kws={"color":"red","linewidth":3}, hist=False, kde=True)

    df_nat[df_nat["Kind"]==tkind].plot.hist(density=True, ax=ax, color="blue", alpha=0.3, label="NatForc")
    if fit_norm: 
        sns.distplot(df_nat[df_nat["Kind"]==tkind]["Trend"].to_numpy(), ax=ax, fit_kws={"color":"blue","linewidth":3}, fit=stats.norm, hist=False, kde=False)
    else:
        sns.distplot(df_nat[df_nat["Kind"]==tkind]["Trend"].to_numpy(), ax=ax, kde_kws={"color":"blue","linewidth":3}, hist=False, kde=True)

    ax.axvline(x=df_obs[df_obs["Kind"]==tkind]["Trend"].to_numpy(), color="black", label="Observed", linewidth=3)

    ax.legend(["AllForc","NatForc","Observed"])

    subtitles= {
       "ALL":"Anomalous Days",
    }

    ax.set_title('Linear Trend in Anomalous Days', fontsize=16)

    ax.set_xlabel('Linear Trend [days per year]', fontsize=14)

    ax.axvline(x=0, linestyle='dashed', color="gray")

    ax.set_ylabel("Probability Density [%]", fontsize=14)


    ax.set_xlim(xranges[tkind][0],xranges[tkind][1])

    plt.yticks(fontsize=14)
    plt.xticks(fontsize=14)

    xi = df_nat[df_nat["Kind"]==tkind]["Trend"].to_numpy()
    yi = df_all[df_all["Kind"]==tkind]["Trend"].to_numpy()

   #percentile rank
    oi  = df_obs[df_obs["Kind"]==tkind]["Trend"].to_numpy()
    if fit_norm:
        mu, std = norm.fit(yi)
        p = np.random.normal(mu,std, 1000)
        obs_all_rank = stats.percentileofscore(p, oi)
    else:
        obs_all_rank = stats.percentileofscore(yi, oi)
    plt.text(0.02, 0.95, "Percentile rank of observations based on AllForc: %.1f%%" % (obs_all_rank), transform=ax.transAxes)
    if fit_norm:
        mu, std = norm.fit(xi)
        p = np.random.normal(mu,std, 1000)
        obs_nat_rank = stats.percentileofscore(p, oi)
    else:
        obs_nat_rank = stats.percentileofscore(xi, oi)
    plt.text(0.02, 0.90, "Percentile rank of observations based on NatForc: %.1f%%" % (obs_nat_rank), transform=ax.transAxes)

    plt.tight_layout()
    plt.savefig("./Fig5.png")
